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Abstract. We propose a description for transient penetration simulations of miscible and immiscible fluid 
mixtures into anisotropic porous media, using the lattice Boltzmann (LB) method. Our model incorporates 
hydrodynamic flow, diffusion, surface tension, and the possibility for global and local viscosity variations to 
consider various types of hardening fluids. The miscible mixture consists of two fluids, one governed by the 
hydrodynamic equations and one by diffusion equations. We validate our model on standard problems like 
Poiseuille flow, the collision of a drop with an impermeable, hydrophobic interface and the deformation of 
the fluid due to surface tension forces. To demonstrate the applicability to complex geometries, we simulate 
the invasion process of mixtures into wood spruce samples. 
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1 Introduction 

Fluid invasion and flow in porous media are ubiquitous 
phenomena in nature and technology. In principle, studies 
on porous media can be classified by three length scales: 
the pore scale, the representative volume element (RVE) 
scale, and the domain scale [T]. Studies on the pore scale 
directly model the pore space geometry and the fluid hy- 
drodynamics [2] . The intermediate sized RVE is the mini- 
mum volume required to characterize the flow through the 
medium. Most systems are studied on the RVE scale, using 
semi-empirical theories. This is necessary due to complex 
pore space geometries. Typical examples are the Darcy, 
the Brinkman-extended Darcy, and the Forchheimer - ex- 
tended Darcy models p]. The flow in porous media can 
also be modeled via generalized Navier-Stokes equations 
[31415] ■ including all necessary terms in one momentum 
conservation equation. Unfortunately in the majority of 
cases, the analytical solution becomes very difficult and is 
even complicated for anisotropic, structured porous media 
like wood. Therefore, one depends on numerical methods 
in this field. The fiuid itself can be quite complicated as 
well. It can be a miscible or immiscible mixture of differ- 
ent species, its viscosity can change locally due to inter- 
nal chemical reactions like hardening and can depend on 
shear velocities in the case of Non-Newtonian fluids. Ex- 
amples are adhesive penetration at wood junctures or the 
impregnation or sealing of concrete for protective reasons. 
A practical model approach for these phenomena is hydro- 
dynamic dispersion [5] that treats one fluid, e.g. a solvent, 
by the hydrodynamic theory and the others, e.g. solutes. 



by diffusion equations. For complicated anisotropic porous 
media and fluids, coupling the diffusion with the gener- 
alized Navier-Stokes equations makes analytical solutions 
quite difficult. However, numerical solutions can be found. 
The motivation for the present work is adhesive penetra- 
tion in wood. 

Wood constructions depend on good adhesive bond- 
ings, that are the result of a gluing process, where hard- 
ening liquids are pressed into a porous, anisotropic mi- 
cro structure. Depending on the penetration and adhesive 
hardening characteristics, preferrable interface morpholo- 
gies can be obtained [7] . The penetration is dominated by 
the strong anisotropy of the permeability tensor and the 
viscosity evolution of the adhesive. Due to the difficulty of 
the problem, penetration simulations with classical con- 
tinuum mechanics are hardly feasible. 

Over the last two decades, lattice Boltzmann (LB) 
methods developed to be an alternative to the simula- 
tion of partial differential equations (PDEs). Originally 
LB methods were developed as discrete realizations of 
kinetic models for fluids [819] . Extensions allow for the 
simulation of diffusion [10], waves |ll|12j . magnetohydro- 
dynamics ^13il4il5] . quantum mechanics [16] . and multi- 
phase flows |17I18I19I20I21I22| . Tolke et. al. made a first 
study on binary flows in porous media |23j . The authors 
consider two immiscible fluids that evolve inside a pore 
scale system by the Navier-Stokes equations with interac- 
tion terms. LB models for incompressible flows through 
porous media on the RVE scale were first addressed by 
Quo et. al. 24J, who proposed a 2D model including a non- 
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linear Forchheimer term, that can be used for isotropic, 
heterogeneous porous media. 

We introduce a 3D LB model that recovers the general- 
ized Navier-Stokes equations in anisotropic porous media, 
extended by the diffusion equation for the solute concen- 
tration in fluid mixtures. Consequently viscosity changes 
appear, that can be superimposed by local hardening. Ad- 
ditionally we incorporate the free surface technique ,25] to 
allow for immiscible fluids. In this paper, Sec. [2] first de- 
scribes some LB fundamentals and the numerical model 
details for fluid and diffusive particles. Sec.[3]shows a brief 
description of the free surface technique while Sec. H] ad- 
dresses the implementation of the algorithm and Sec. [5] 
demonstrates the validity of our model and its applica- 
tion to diverse physical systems. 



2 3D Lattice- Boltzmann Model for 
Anisotropic Porous Media 

We begin with the generalized Navier-Stokes equations for 
low Reynolds numbers with mass conservation 

^+V-ipV) = , (1) 
and momentum conservation, 

^ + (V-V)(X) = -^ + ..V^V + .G (2) 

p ' 

with the fluid density p and viscosity v, the volume-averaged 
velocity and pressure V and P. e denotes the porosity 
of the medium, K the permeability tensor that can be 
anisotropic, and G represents an external force field. The 
quantity is named the effective viscosity which is not 
expected to be the same as the viscosity of the fluid due 
to the effects of tortuosity and the dispersion of viscous 
diffusion flux [I]. The surface tension tensor S can be ex- 
pressed by 

^'^1 V F^F, F^Fy -F^-F^J 

with the color gradient F=Vp, a parameter A that flxes 
the strength of the surface tension, and the characteristic 
relaxation time r, that is flxed by the effective viscosity 
i/e=\{T — i) [53]. Most of these quantities are related to 
each other like the pressure P=pc^g/e depending on the 
speed of sound of the fluid c^. In the following we demon- 
strate how Navier-Stokes equations are reproduced by LB 
methods, and how extensions like external forces, Darcy 
flow, and surface tension can be added to the framework. 
Diffusion however is treated differently and explained sep- 
arately. 




Fig. 1. Cubic lattice D3Q19 discretization of the fluid velocity. 
The arrows represent the velocity vectors vf , where p indicates 
the plane of location. Note that by including the rest vector 
located in the center of the cell, totally 19 directions are de- 
fined. 

2.1 Model Description for the Generalized 
Navier-Stokes Equations 

The basic principle of LB is that Navier-Stokes equations 
are not solved directly on a grid of cells, but via a strongly 
simplifled particle micro dynamics with discrete Boltz- 
mann equations and collision rules such as Bhatnagar- 
Gross-Krook (BGK) 27J . The two ingredients particle stream- 
ing and collision are calculated along a certain number of 
fixed orientations (19 in our case as illustrated in FigU]) in 
each cubic cell with lattice constant 5x=5t. The velocity 
vectors are denoted by v^, where z=1...6, indicates their 
orientation and p=0, 1, 2 their reference plane (see FiglT]). 
Each velocity vector comes along with one assigned dis- 
tribution function ff . The density p in a cell is the sum 
over its distribution functions 

P-ttff • (4) 
We obtain continuum properties for V and F via 

/'V-EE/M ' F = EE/^(^+<K ■ 

i—l p—0 i—1 p—0 

In order to include the influence of the porous media, of 
external forces, and of surface tension on the momentum 
equation (Eq.[2]), we need to include additional force terms 
in the LB model. Following the work of Guo et.al. [55], we 
add two force terms T[ and Tq to each of our 19 directions, 
including the rest vector in the cell's center: 

/f (X + vf, i + 1) - /f (x, t) = I2f (X, t) + T[, (6) 
/o(x,t+l)-/o(x,f) = l2o(x,t)+To . 

The second ingredient of LB are the BGK collision terms 

nf and Ho [m 

= -i(/f(x,i)-/r^(x,i)) , (7) 
no --i(/o(x,t)-/o"^(x,i)) . 
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The subscript '^'^ denotes the equihbrium states for the 
distribution functions, which we must find such that the 
model reproduces correctly the hydrodynamic equations, 
r denotes again the relaxation time from Eq. [3l The force 
terms T[ and Tq are given by [28j : 

T[^(^- ^) . F-*) + f K • V')(vr • F-*)(8) 



,3 (^Y' • F"^^*) 
To ( 1 - J- ) ( -^V • F^^* 



with the external force 



(9) 



in our case. Due to the external forces, the velocity V 
needs to be corrected to obtain V' by 



V' = V + -W^^* 



(10) 



Unfortunately F''^' depends on V'. Then we insert Eq. [9] 
into Eq. and obtain 



V = • ( V + 



(11) 



where we define A=l + ^K^^ with the identity matrix I. 

We finally write down the equilibrium functions for our 
system: 



l + 3(vf .V') + ^(vf .V')2 (12 



_3.v'2 , Act I (vf'F)" 
'2e^ P IFI 



2e 



AerwolFl 



(13) 



with the known weights wq — 1/3, u'i_2.3,4 — 1/36, and 
u;5,6 = 1/1 8 for the D3Q19 cells |29]- The last terms in 
Eqs. I12I13I represent the surface tension. In the contin- 
uum limit, the surface tension tensor S from Eq. [2] is 
reproduced. Note that it can be shown analytically [30] 
that this LB BGK model recovers the generalized Navier- 
Stokes equations Eqs. I1I2I in the isothermal and incom- 
pressible limit with surface tension. 



2.2 Lattice Boltzmann model for solute diffusion 

To consider diffusion in an LB framework, we can either 
define an additional set of velocity vectors and include 
their diffusive terms directly into the equilibrium func- 
tion (Eq. [T2|) , or define additional equilibrium distribution 
functions for the diffusion but work on the same velocity 
vectors as proposed by Hiorth et. al. [31] on a D2Q9 cell 



configuration. In this work, the second approach was cho- 
sen, to simplify a later extension to various species. We 
therefore recover the diffusion equation in porous media 
based on Ref. [31], but extend it to heterogeneous porous 
media. Starting from the diffusion equation 



dC_ 
'dt 



V • (CU) ^ V • [DVC] 



(14) 



with the concentration C, the diffusivity D, and the ve- 
locity of the mixture U, Hiorth et. al. [31i use a D2Q9 
cell configuration and define a distribution function hi for 
each velocity vector. The concentration is calculated via 
C—Y^^^^ hi and the equilibrium equations are defined by 



l + 3(vf -U) 



WqC 



(15) 

with the weights wo=4/9, wi, 2, 3, 4=1/9, and W5. 6,7, 8=1/36 
|31) . To extend to 3D and to include the porous medium, 
we propose modified equilibrium distribution functions on 
a D3Q19 cell configuration (see Fig. [Ij 



l+3(vf-V* 



/iq (x, t) = wqcC 



(16) 

where Y*=V'+D^. With the correct weights uiQ^i for the 
D3Q19 cell, we evolve Eq. [16] according to the Boltzmann 
equation with the BGK collision term Sl^^ (x, t) [27] : 

;.f(x + v„t+l)-/.f(x,t) = /2^,(x,t) , (17) 
/2,,.(x,t) = -^(/.f(x,t)-/.rV,i)) ■ 

Now Tjj defines the relaxation time for the diffusion model 
with the diffusivity, I?=^(t£) — ^)- We can prove via a 
Chapmann-Enskog expansion that the model reproduces 



dt 



V • (eCV) = V ■ [eDVC] 



(18) 



in the continuum limit. 

The full model can now reproduce the generalized Navier- 
Stokes equations in the continuum limit for anisotropic 
porous media, including surface tension and the diffusion- 
advection of the solute. To complete the description of the 
model, we will briefly summarize the free surface technique 
[25] used to describe the movement of fluid with a fluid-gas 
interface. 



3 The Free Surface Technique 

If a fluid penetrates an unsaturated porous medium, free 
surfaces form. To include free surfaces in our LB scheme 
a free surface technique [25] is necessary. If a cell is only 
filled by gas with a negligible density, Eqs. [T^ and [T51 lead 
to vanishing equilibrium functions and the simulation be- 
comes unstable. This problem can be solved by classifying 
fluid cells into three types of cells: the liquid cells, to- 
tally filled by liquid fluid, the empty cells entirely filled by 
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gas, and the interface cells that contain both fluids: liq- 
uid and gas. The purely liquid filled cell were discussed in 
the previous section, while gas filled cells can be excluded 
from the further calculation due to the high density ratio 
between fluid and gas. However the interface cells need 
further attention. 

Interface cells consist partially of liquid and of gas, 
determined by the fluid fraction A(x, i) = rn(x, f)//9(x, t) 
with the mass m of liquid in the cell located in x at time t. 
Using a normalized cell volume, the mass for the fluid cell 
is equal to the density of the liquid, and for gas cells equal 
to zero. Analogous to Eq. [H we obtain the macroscopic 
mass m(x, t) by summing the mass functions over all 
directions 

m(x,i) = ^^mf . (19) 
The mass functions evolve according to 

mf(x,i + l)=mr(x,t)+^^zimf(x,i+l) , (20) 
with 

Am^i., t + l) = /;(x + vf , ,^M^i±<A±M^^,,) 

_ ,, A(x + vf,t) + A(x,t) 
Ji l^i 2 ' 

where the index j stands for the vector with opposite di- 
rection to the vector with index i. The interface cells also 
have associated distribution functions that evolve via the 
discrete Boltzmann Eq. [H] identical to those in the fluid 
cells. These distribution functions are used to calculate 
the density and the updated fluid fraction A. 

Nevertheless, we still obtain zero distribution functions 
in the interface cells from Eq. [6l due to empty neighboring 
cells. To avoid this problem, we calculate the evolution of 
distribution functions of interface cells originating from 
empty neighboring cells with the modified function 

/;(x,t+l) = /r'^(x,t) + /f^V,i)-/f(x,<) . (22) 

During the system evolution, the interface moves, leading 
to changes in the fluid fraction A. If A becomes unity, the 
cell changes to a liquid cell, while A ^ leads to gas or 
empty cells. To complete the picture all cells neighboring 
new liquid cells are set to interface cells. To summarize, 
by classifying cell types and corresponding distribution 
functions, depending on fluid type, we can capture the 
evolution of free surfaces. A more detailed description can 
be found in Ref. [25j. Finally this technique extends our 
model to simulations of fluid penetration with gas/liquid 
interface. Before validation examples are shown, we out- 
line the implementation of the algorithm. 

4 Algorithm Implementation 

To obtain a transparent, portable simulation code, we im- 
plemented it object oriented in C++. The general frame- 
work of a simulation is given in Fig. [2] First we set all 



[ Set initial fields: p, V, C. (Eqs. HIS)) 

Define fluid, interface and empty cells (free sur- 
face technique 25 ) and set the mass of each cell. 

Calculate: /f and h^''^^ from fields, 
/f^/f^'^ andfef=fef^+ (Eqs.imaiSl) 

^ i 

[ t=t+i. 

I 

\ 

Stream from all cells to neighbor- 
ing cells: /f(x+vf,t+l)=/f(x,t) and 
ftf(x+vf,t+l)=/inx,t). (Eqs.iHTl) 

Stream of the mass in all cells to neighboring 
cells (free surface technique) 25 . (Eqs|20) 

I 

• \ 

Calculate fields taking moments 
of /f and /if in aU cells. (Eqs. gEJ 

^ Calculate the color gradient F. (Eq. [5| ] 
I 

Apply collision in all cells: 

/f(x,^+l)=/f(x,^)-^?f(x,^)+Tf and 

/if (x, t+l)=/if (x, f)-^?f>,(x, t). (Eqs. 

Redefine fiuid, interface and empty 

cells (free surface technique) accord- 

ing to the mass in each cell |25| . 

Fig. 2. Computational procedure of the LB simulation of 
fiow of mixtures through anisotropic porous media including 
free surface and tension effects. 



initial field values like the density p, the velocity V, ex- 
ternal force field F^^* (Eqs. I4l5|) . and the solvent concen- 
tration C inside the mixture. To completely define the 
system, the porosity field, the permeability fields, and 
relaxation times are set. Also the initial saturation and 
boundary conditions like periodicity, pressures, flow rate, 
or impermeability are imposed. We can define three types 
of cells: fluid, interface and empty, depending on the prob- 
lem. With these quantities set, we can start calculating the 
distribution functions /f , from the equilibrium distri- 
bution functions ff''"^ and /if^"' fEas. 112113116)1 . Now, the 
time step can be incremented and system is ready to evolve 
(see Fig.O. 

To follow the evolution of the system, first solve the 
stream to all neighboring cells, stream their mass and cal- 
culate the macroscopic variables. Then the color gradient 
F is calculated. Finally, we have to update the type of 
cells and increment the time to prepare the system for the 
next time step. These steps are looped over all cells in 
the system. We stop the simulation, when a stable state 
is reached. 
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5 Model Validation and Calculation of the 
Flow Profile of Adhesives in Wood 

To validate our LB model implementation, we simulate 
systems with known analytical solutions that especially 
address each term that appears in the generalized Navier- 
Stokes equations. In particular we simulate Poiseuille flow 
through isotropic and anisotropic porous media, droplet 
formation, and fluid surface smoothening. Finally, we show 
an application for the penetration of hardening mixtures 
into an anisotropic, heterogeneous porous medium, demon- 
strating the capability of the model to represent adhesive 
penetration in complex materials like wood. 

5.1 Single Fluid in Homogeneous, Isotropic Porous 
Media 

To show the capability of the model to reproduce the 
Darcy law in a simple case, we simulate a generalized 
Poiseuille flow driven by a constant force in an isotropic, 
homogeneous medium between two infinite plates at dis- 
tance L (see FiglS]) . For this case, we need neither the sur- 
face tension tensor nor the free surface technique. There- 
fore just the force terms from the Darcy law, body force 
and the viscosity terms on the right hand of Eq. [2] are 
active. We impose periodic boundary conditions on all ve- 
locity vectors and null velocity at y=0, L of our 3D sys- 
tem. We only consider movements in the a;-dircction, so 
the velocity has form V=(T4, 0, 0) with boundary values of 
Vx{x, 0)=Vx{x, L)=0. The analytical solution of this prob- 
lem is well known [I] by a Brinkman-extended Darcy equa- 
tion, and has the form, 



where r= 




We set up a system of 3x100x1 cells with porosity 
e = 0.3, relaxation time t = 0.6, density p—1.0, and the 
body force G=(0.1,0, 0). The permeability K is fixed by 
the Darcy number Da=K/L^ to lO^^^, SxlO^^^, and 10~* 
respectively. We start with a flat profile and let the sys- 
tem evolve. When the velocity profile reaches a steady 
state, the simulation is stopped. Comparing the solution 
to Eq.[23]we find excellent agreement with the theory (see 
Fig.©. 

5.2 Anisotropic Permeability 



° Da =1x10" 
o Da = 5x10 




20 40 60 80 100 

y (cells) 



Fig. 3. Velocity profile for a generalized Poiseuille flow in 
a homogeneous, isotropic porous medium for various Darcy 
numbers. The generalized Poiseuille flow is obtained when the 
fluid moves through a channel between two walls driven by a 
constant body force through the porous medium (shaded zone). 

X (cells) 

10 20 30 40 50 60 70 
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□ V^(x)/V^___^,^ 






Eq. 20 
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Fig. 4. Velocity profiles for an anisotropic porous medium 
with diagonal permeability tensor. Squares denote the ve- 
locity profile of the x component and diamonds represent 
the simulation with flow in y-direction {vx,max ~ 3 • 10""^, 

^y,max — 3 " 10 ). 



, , GKf^ cosh[r(y - L/2)]\ 

Ky- 1 ^ Z T /9 ' 23 0.6 

V \ coshrL/2J / 



In order to validate the model for anisotropic porous me- 
dia, we study Poiseuille flow like in Sec. 15. li but with a 
diagonal permeability tensor K, that was fixed to K = 
10-34, + 10-5^^^. 

First, we let the fiuid move inside the porous medium 
in x-direction like before with T4(y = 0, L)=0 at the bot- 
tom and top. In a second simulation the fluid moves in z- 
direction with Vy{x = 0, L)—0 at the vertical boundaries. 



The velocity proflles can be calculated again by Eq. [23l In 
the simulation we used an array of 70 x 70 x 1 cells with the 
same porosity, relaxation time and density as previously. 
However the body force is G=(0.1,0, 0) for the first case 
and G=(0, 0.1, 0) for the second one. We run the simula- 
tion until the system reaches a steady state and measure 
the respective velocity profiles (see Fig.Hl). The simulation 
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2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 

Time (clicks) 

Fig. 5. Evolution of the total mass of a falling droplet that 
crashes against an impermeable, hydrophobic wall. Inset: Snap- 
shots of the droplet at initial time (a) and after 20000 time 
steps (b). 

results show that the model reproduces the Darcy flow in 
anisotropic porous media accurately. 

5.3 Free Surface Model 

As condition for a vaUd implementation of the free surface 
technique, we have to assure, that our model conserves 
mass. This can be tested by simulating a freely falling 
droplet followed by its collision against an impermeable, 
hydrophobic wall. 

The system consists of a grid of 140x140x1 cells. Ini- 
tially a droplet with a radius of 8 cells is placed in the 
upper part of the system and an impermeable surface 
avoids its penetration through the bottom of the simu- 
lation zone at y=0. The system parameters are fluid den- 
sity p = 1.0, relaxation time t=0.7, external gravitational 
force G=l • 10~^ • (0, 1, 0). The impermeable wall is mod- 
elled by setting the macroscopic velocity in the respective 
cells to zero. The force terms from the Darcy law and 
surface tension are neglected here. 

In Fig. [5] we plot the droplet mass as function of time. 
We can see that the mass is conserved with a small er- 
ror of less than 1%. Note also that this error is within 
the range of common errors that go along with LB simu- 
lations. After approximately 10600 time steps we flnd an 
abrupt peak in the total mass due to the collision with the 
wall. Note that due to numerical reasons, for a short time 
mass is transferred to the boundary cells that represent 
the impermeable wall, but is returned in the next time 
steps, so that the total mass recovers the value previous 
to the collision. The initial and flnal conflgurations after 
20000 time steps arc shown in the figure's inset. 

5.4 Introduction of the Surface Tension 

Mass conservation alone is not sufficient to describe free 
surfaces, rather surface tension needs to be considered. A 




Fig. 6. Surface smoothing due to surface tension from the 
initial shape (a) to the final shape after 10000 time steps (b). 
Intensities represent the quantity of fluid in each cell where 
black means empty cells and white completely filled ones. 

correct implementation can be for example tested by sim- 
ulating the deformation of the fluid from a rectangular to 
a circular shape. In this spirit we implemented a simula- 
tion using an array of 100x100x1 cells, with an initially 
flUcd, square fluid zone of size 20x20 cells (see Fig. 15]). The 
density and relaxation time are identical to the droplet 
simulation, and a surface tension parameter ^4=0. 001 was 
chosen (see Eq. [3]). The force terms from the Darcy law 
and external forces like gravity were switched off for this 
simulation. The system evolved until a steady state was 
reached after 10000 time steps. The final shape of the fluid 
surface smoothened due to surface tension is seen in Fig. [6] 

Summarizing, we showed that the model can repro- 
duce the theory of single fluids in the case of anisotropic 
porous media with surface tension effects and external 
forces. Note that in our case it is only possible to ap- 
ply the free surface technique to liquids with low Reynold 
numbers. 

5.5 Penetration of Adhesives in Wood 

To demonstrate the applicability of our model, we describe 
the simulation of the penetration of adhesives into wood 
which is an anisotropic porous medium. The adhesives 
can be treated like two fluids: a polymer (solute) and a 
solvent, e.g. water. The polymer dynamics is governed by 
the diffusion equation and the solvent by the generalized 
Navier-Stokes equations. To consider hardening effects in 
the adhesive, we have to consider a local change in time 
of the viscosity due to the change of concentration. An 
additional difficulty is the description of the anisotropic 
porous medium wood. We need to obtain the permeability 
tensor of wood and to implement a law for the dependency 
of the viscosity with the concentration and time. 

Wood properties: In wood, the orthotropic permeability 
tensor and density is a local property, depending mainly 
on year ring geometry, orientation, and position inside the 
year ring. Following Refs. |32I33| an empirical expression 
of the wood density of the type 

p(a;) = po[l + 2aa::''exp(-ca;)] , (24) 
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with the parameters a, 6, c, the position x, and the mini- 
mum density po can be made with parameters depending 
on the type of wood. For spruce we can take values of 
Po=310kg/m3, a=^, 6=0.8, and c=^, where Z\=2.9mm 
is the year ring width [32] ■ To represent several paral- 
lel year rings, we use a periodic function of the form of 
Eq. [21] with periodic A. Since wood consists of cell walls 
and lumen, we can estimate the porosity from the density 

by 

e{x) = l-f^ , (25) 

Pmax 

where Pmax is the density of cell walls. Permeability and 
porosity are related in the scalar case by [1] 



K = 



n{l 



Z?2 

e)2 P 



(26) 



with the characteristic pore size Dp and a constant n that 
parametrizes the microscopic geometry of the material 
taken as n=150 in our case. In wood, the pore size Dp 
must also depend on the wood density, and we propose as 
a first approximation the relation 



Dpix) 



Co 1 



p{x) 



(27) 



with a constant C2 , that can consider the shape of the mi- 
croscopic pores [31]. Inserting Eqs. [25] and [23 in Eq.l^we 
obtain the macroscopic permeability of a sample as a func- 
tion of the year ring coordinate x. C2 can now be obtained 
by integrating Eq. [211 over one year ring and comparing 
it to literature values for the macroscopic permeability 
|35I36| of spruce in radial, longitudinal, and tangential di- 
rection. We finally calculate the required local permeabil- 
ity tensor in cylindrical coordinates (r, z, 9), as 



K = 



p{x) 



e{x)- 



8.2 
51470 
0.82 



(28) 



Hardening Process and Viscosity Model: To determine 
the viscosity evolution during the hardening process, we 
have to take into account two main factors: the chemical 
reactions that take place between the adhesive molecules 
and the solute concentration. The first factor causes a vis- 
cosity increment with time, while due to the second one, 
viscosity changes with local solvent concentrations. In or- 
der to describe these mechanisms, we choose a viscosity 
model given by 



i/(C, t) ^i'g[l + aexp(at)] exp(/?[l - C]) 



(29) 



where Vg^ a, a, and (3 are parameters that can be obtained 
from experimental hardening curves of adhesives with di- 
verse solvent concentration. 
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Fig. 7. Fluid profile during the penetration into the wood 
structure. Intensities are identical to Fig. [6] The wood sample 
is located in the region j/<60 with a year ring inclination of 
45°. 



Simulation of Adhesive Penetration into Wood: Em- 
bedding the permeability tensor field and the hardening 
function, our model is ready to simulate the penetration 
of adhesives into the porous wood mesostructure. We use 
an array of 90x90x1 cells and locate the wood sample 
in the zone ?/<60. The year rings are inclined by 45° to 
the horizontal plane. We take the fluid density p — 1.0, 
the parameters for the porosity, according to Eq. [211 with 
a=0.62, 5=0.2, and c=— 0.2, corresponding to spruce wood 
[32] . and a year ring distance A of 30 cells. The viscosity 
model uses the parameters Ug=Q.\^, a=2.05 x 10~^, and 
/3=10~'^, obtained from viscosimetric measurements with 
urea formaldehyde adhesive [37j . Additionally, we applied 
a constant external pressure of Pea;t=0.0441 that is equiv- 
alent to 0.9mN/mm^ in IS units. 

The simulation was stopped when the adhesive could 
not penetrate further due to hardening. Fig. [7| shows the 
penetration of the fluid into the wood sample. We can 
see that the model can simulate complex materials like 
wood. In some zones the adhesive penetrates with more 
speed due to the higher local porosity and the direction of 
the movement is correlated with the principal axis of the 
permeability tensor. 



Discussion and Conclusions 

We introduced a new lattice Boltzmann model to simulate 
the dynamics of the flow of mixtures in anisotropic, het- 
erogeneous porous media. Our 3D model can be applied 
to many problems in material science and engineering. It 
includes a free surface technique to simulate the invaded 
fluid profile inside the material structure and a surface 
tension term to control the inerface shape dynamics and 
the forces of cohesion inside the fluid. 

The accuracy of the model was tested for a set of sim- 
ple cases, like generalized Poiseuille flow in isotropic and 
anisotropic porous media, droplet formation and surface 
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smoothing. Wc found excellent agreement with the Darcy 
law for Poiseuille flow. The surface free technique vali- 
dated in the simulation of a freely falling droplet and its 
collision against an impermeable and hydrophobic surface 
showed negligible errors. The surface tension effects were 
tested by simulating the smoothing of the fluid profile. To 
demonstrate the applicability of the model to real cases, 
we implemented a simulation of the penetration of adhe- 
sives into a sample of spruce wood. We showed that the 
model can reproduce the dynamics of hardening fluids in 
complex materials even when the media are anisotropic 
and heterogeneous. 

The actual model has all the advantages known for the 
Lattice Boltzmann method like the minimal use compu- 
tational resources, fast algorithms that allow for real time 
simulations, and the ability for efflcient parallelization. 
Additionally extensions of the model are rather straight- 
forward. We hope that the proposed simulation approach 
will be useful in order to model complex geometries in het- 
erogeneous and anisotropic porous media even for compli- 
cated fluids like mixtures of two components. 

The authors arc grateful for the financial support of the Swiss 
National Science Foundation (SNF) under grant no. 116052. 
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